Van der Waals Energies in Density Functional Theory 



On 



O 

CO 



> 

00 
CN 

m 
o 



c 

o 
o 



13 



Walter Kohn 1 , Yigal Meir 2 ' 3 and Dmitrii E. Makarov 4 

1 Department of Physics, University of California at Santa Barbara, Santa Barbara, CA 93106 
2 Physics Department, Ben Gurion University, Beer Sheva 84105, ISRAEL 
Institute of Theoretical Physics, University of California at Santa Barbara, Santa Barbara, CA 93106 
4 Department of Chemistry, University of California at Santa Barbara, Santa Barbara, CA 93106 

In principle, density functional theory yields the correct ground-state densities and energies of 
electronic systems under the action of a static external potential. However, traditional approxi- 
mations fail to include Van der Waals energies between separated systems. This paper proposes a 
practical procedure for remedying this difficulty. Our method allows seamless calculations between 
small and large inter-system distances. The asymptotic H-He and He-He interactions are calculated 
as a first illustration, with very accurate results. 



Density functional theory (DFT) 0] has become a use- 
ful tool for calculating ground-state energies and density 
distributions of atoms, molecules and solids, particularly 
of systems consisting of many atoms j|. The simplest 
approximation for practical purposes is the local-density 
approximation (LDA) H , based on the properties of the 
uniform electron gas. The so-called generalized gradient 
approximations (GGA) are important refinements of the 
LDA. 

In principle DFT yields the exact ground state energy, 
including long-range Van der Waals (VdW) energies, very 
important in organic chemistry and elsewhere. However 
the commonly used LDA and GGA, designed for non- 
uniform electron gases, fail to capture the essence of VdW 
energies. The latter reflect correlated motions of elec- 
trons due to the Coulomb interactions between distant, 
even non-overlapping atoms, molecules and solids. Thus 
a new strategy is needed. 

Here we propose a first-principles approach, which con- 
tains the following essential ingredients: 1. The density 
distribution, n(r), is approximated by the LDA or GGA. 
2. The Coulomb interaction is divided into short and 
long range parts, of which only the latter contributes to 
VdW energies. 3. The contribution of the long-range 
interactions to the energy is expressed by the so-called 
adiabatic connection formula (see Eq.(^) below). 4. This 
expression is transformed into the time-domain, avoid- 
ing the need to solve a self-consistent equation for the 
density-density response function. As an illustration we 
calculate the asymptotic VdW interaction between two 
Helium atoms and between Hydrogen and Helium atoms, 
with excellent results. The method allows seamless cal- 
culation of the interaction of two subsystems, e.g., an 
atom and a surface, from small to large separations. Our 
work was carried out independently and differs substan- 
tially from recently published work by Anderssen et al 
H and by Hult et al. which depend critically on a 
fitting parameter. 

Since the VdW energies are due to the long range of 
the electron-electron interaction, U(r) — 1/r, we sepa- 
rate this interaction, as a preliminary step, into short 
and long range parts, 



U{r) = U sr {r) + U tr {r). 



(1) 



electron-electron distance, so that the effect of Ui r on the 
total energy is small. The calculated total energy is, in 
principle, independent of the choice of K, in practice — 
with appropriate approximations — nearly so. 

We now write the Hamiltonian as a function of a cou- 
pling constant, A, that "turns on" Ui r , such that the 
physical Hamiltonian operator corresponds to A = 1: 



H(X) = T + V x + U sr + XU ir , 



0<A<1 



(2) 



where T is the kinetic energy, and the external potential 
V\ is chosen such that the ground state density ri\(r) of 
TC(X) equals the exact physical density n\ = i(r) for all A 
[Q. Note that for A = 0, the interaction is entirely short 
range and that for A = 1, V\ = \ — V e xt- We denote the 
ground state energy of H(X) by E(X). Then the ground- 
state energy of the physical system, E = E(l), is given 
by 

E = E(0) + fdr [V cxt (r) - V (r)} n(r) 



(3) 



+ ~ / .InlrT ,,(,■- /•') 



(h(r)h(r ))x dX — n(r)5{r — r ) 



where h is the density operator; Vo(r), defined above, 
eventually drops out; see Eq.(||). 
From DFT, E(0) is given by 

E(0) = T s [n(r)] + [ drV (r)n(r) 



+ - J drdr' U sr (r - r')n(r)n(r') + E s x r c [n(r)}, (4) 

where T s [n(r)], is the non-interacting kinetic energy func- 
tional and E^ r c [n(r)] is the exchange correlation energy 
of an electron gas with density n(r) and the short range 
interaction U sr . 

Substituting (|J) in (Q) wc find, after simple manipula- 
tions, the exact result (independent of the form of Ui r ): 

E = T s [n{r)] + [ dr V ext (r) n(r) 



For example, we can choose U sr (r) = e Kr /r, with 
k _1 chosen somewhat larger than a typical intra-atomic 



+ — J drdr'U(r — r') n(r)n(r') 
+ EZ\n(r)] - U ir (Q)N + E pol [n(r)} (5) 
where N is the number of electrons, and 



1 



E poi = 2 / drdr' Ui r {r - r') 

xf dX((h{r)-n(r))(h(r')-n(r')))x . (6) 
Jo 

£"p i includes the long-range polarization energies. 

To calculate the first four terms in Eq. (||) we use tradi- 
tional methods. Experience j^j has shown that the den- 
sity n(r) may be calculated to a very good approximation 
by the LDA with the full U(r). Such a calculation also 
automatically yields an approximate T s [n.(r)], 



N . 

T sM = ^2 \ v '<s( r ) n ( r ) dr, 



(7) 



where vks is the Kohn-Sham (KS) effective poten- 
tial that reproduces n(r), and the tj are the single- 
particle energies available from the LDA calculation. For 
E^, r c [n(r)) there exist unpublished excellent results in the 
LDA g. 

To evaluate E po \ we use the appropriate exact connec- 
tion formula ||, 

r rl r<x> j 

Epoi =-J drdr'U er (r - r') J dX J ^Imx(r, r'; u, A), 

(8) 

where the linear-response susceptibility \ is defined, as 
usual, as follows. Let aVi(r, uj)e~ luJt be a small perturb- 
ing potential, acting on the ground state of 7i(A), and 
producing a density response an\{r 1 w)e~ Mt , where a is 
infitesimal. Then x is defined by 



ni(r,oj)= dr' x(r,r';ui,X)V 1 (r',uj) 



(9) 



For A = 1, methods to evaluate x( r > r 'i ^, A) have been 
discussed in the past f|||, and can be formally carried 
over to A < 1. x is the solution of the integral equation: 



X (r, r'; uj, A) = Xi<s{r, r';uj)+ J dr" ' dr"' 'x KS (r, r";uj) 

x \U{r" - r>") + f xc (r", r'";oj, A)] X (r>" , r>; to, A), (10) 

where Xks{t, r';u>) is the response function of the cor- 
responding non-interacting Kohn-Sham system, and f xc 
describes exchange and correlation effects (see Eq.(6) of 
Ref. |). 

However, except for systems of very high symmetry, 
such as spherical atoms, the self-consistent solution of 
( fi(i| ) is computationally very forbidding. Here we pro- 
pose an equivalent but much less cumbersome proce- 
dure, which avoids the solution of a self-consistent in- 
tegral equation for each value of A. We note that 
X{r, r';w,X) is the Fourier transform, x( r i r '; w i A) — 
J dtx(r,r';t, X)e luJt , of the time-dependent response 
function, x(r, r';*, A), defined as follows: 

m(r,f,A)= / dr'dt'x(r,r';t-t',X)V 1 (r',t') , (11) 



where Vi(r',t) and n%(r, t,A) are, respectively, external 
perturbing potential and density response. Eq.(||) can be 
rewritten as 

If f 1 f°° dt 

E pol = -—J drdr'U tr (r-r') J dX J - x (r, r'; t, A). 

(12) 

Following Gross and Kohn || , we can replace the density 
response of the physical system to the external perturb- 
ing potential by the response of the (A— independent) 
KS system to an effective potential, 

m(r,t,X) = J dr'dt' XK s(r,r';t-t',X)V 1 eff (r',t',X) , 

(13) 

where 

V l ff {r', t', A) = V x {r', t') + V hxc {r' , t', A) 

+ / dr" [U sr (r' - r") + XU ir {r' - r")] m(r",t, A) . (14) 



V ltXC (r',t',X) is defined by Eqs.(|lj) and (J1J). Thus 
x(r,r';t, X) is the density response of the non- 
interacting KS system at time t, to the V* [r' , t' , A) 
at < t' < t, induced by an external potential 
Vi(r",t') = 5(r" - r)S{t). 

To complete this procedure we need a practical ap- 
proximation for Vi. xc , for which, following Ref, ||, we 
set 



V^ xc {r\t',X) = 



dV xc (n, A) 



On 



m (/,*', A), (15) 



n (r') 



where no is the unperturbed density and V xc is the static 
exchange-correlation potential in the LDA. Here, in ad- 
dition to the usual approximation of the LDA, the fre- 
quency dependence (or retardation) of V xc is neglected. 

The evaluation of X now requires the calculation of 
the evolution of the non-interacting KS system under 
the action of V± (r' , t', A). At this point it is conve- 
nient to change from the coordinate representation to an 
orthonormal basis, and write generically for any 

F(r) and G{r,r') 



F(r) = J2 F mfm(r) 

m— 1 

oo 



(16) 



m,m / =l 

Thus Eq.([ll|) becomes 



A)= / dt'J2xK S ,mm'(t-t';X)V 1 e "(t'^)- 
Jo i 



(17) 



The following steps need to be carried out for each value 
of A: 1. At time t — 0~ the KS system is given by 
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the determinant (7V) _1 Det \(f>i 4>2 ■ ■ ■ 4>n\ of the occupied 
KS orbitals <f>j. 2. At time + , after the action of 
a small external perturbation, V(r) = ctf m (r)5(t) (a 
small), each of the KS orbitals is changed, <f>j(r,t) — > 
(f>j(r,i) — iaf m (r)4>j(r,t). (Effects on the orbitals of 
the finite unperturbed KS Hamiltonian and of the in- 
duced parts of V e M in the infinitesimal interval (0 _ , + ) 
are negligible). 3. For t > + we integrate the time- 
dependent Schrodinger equation for each <f>j in a step- 
wise fashion, evaluating the first-order induced density 
ani m (r,t) = J2j \4 } j{ r ^)\ 2 — n o( r ) at each time step, to 

be able to compute the induced parts of (Eq.(0)), 
which depend on n\ im {r, t). 4. The projection of ni jTn on 
f m > gives Xmm'(t), 



\m m' iP) {fm' •> ^l,m(^)) ■ 

From Eq-IJ), we obtain 

1 f°° dt f 1 > 

Ep i = — — / — dX \ Xmm'(t, X)Ue r 
47T J t J Q ^ 



(18) 



(19) 



In practice, the integration over A is replaced by a finite 
sum. 

As a simple example of this general procedure, we now 
apply it to the calculation of the asymptotic VdW in- 
teraction of a pair of spherically symmetric atoms. We 
denote the atoms by A and B, and their nuclear coor- 
dinates by Ra and Rb (taken to be on the z-axis), and 
write R — \Ra — Rb\- We take R 3> + ag, the 
sum of the atom radii and k ~ (Ra)~ 1/2 , where a is the 
atomic radius. The asymptotic VdW interaction is ob- 
tained from those parts of E po i (Eq.@) in which r and 
r' are in different atoms. Take r to be in A and r' in 
B. Next we write U — U sr + \Uf r . Since k — > when 
R — ► oo, Uir can be treated as a small perturbation, 
giving to first order 



X(r,r';u,X) = A / dr 1 dr 2 U tr (r 1 - r 2 ) 

x XA(r,n;u})xB{r2,r';uj), (20) 

where Xa(Xb) is the response of the isolated atom A (£?). 
The integration over A is now trivial. Lastly, we ex- 
pand Ui r (r — r') in 1/R and obtain the final expression, 
EvdW = —Cq/R , 

3 
3 



C 6 = -Im/ du X z A z (u)x% z (u) 

'0 

X z A z {ti)x z B z {t2) 



dtx 



dto 



(21) 



In the above \ zz i s defined as the z-component of the 
density response to a perturbation in the z-direction, 
X zz = J dr\ dr 2 x( r ii r 2) z\ Z2- The first form is well 
known, the second its Fourier transform into the time 
domain. 

We have calculated the time-dependent response for 
the Helium atom in DFT as follows. We begin with the 



exact V xc {r) Q, which reproduces the exact ground- 
state density no( r ) (known from highly accurate inde- 
pendent calculations), and the corresponding exact KS 
ground-state wavefunction </)q and energy eo- We take 
as perturbation Vi(r,t) — —azS(t). At time t = + the 
wavefunction will be 



(r, t = + ) = 0o (r) - iazcj) (r) 



(22) 



a combination of s- and p-like functions. For t > 0, 
we solve the time-dependent Schrodinger equation for 
4>(r,t) = <po(r,t) + a4>i(r,t), with the initial condition 
Linearizing in a gives the following equation for 



i 9 * 1 ^ = (r, t) + Hi (t)Mr, t) 
<f>i(r, + ) = -ir<j> {r) 



(23) 



e (f>o(r), Ho is the KS unperturbed 
Helium Hamiltonian, 



where (f>o(r,t) — ° leat ' 
, ni(r',t) 



Hi{t) = J dr' 



dV x , 



dn 



»i(r,t), (24) 



n (r) 



and m(r,t) = 4i?e [4>o(r, t)(/>*(r, t)]. V xc was calculated 
using the parameterization of Vosko et al. [^3|. This 
equation was solved by stepwise integration in time. The 
time evolution from t to t + At was carried out using 
the fast Fourier transform method as used in Ref. [ pT[ . 
Since at each instant <j)\ evolves under the action of the 
total effective potential, the resulting response function 
X{t) (and, if desired, the corresponding x( w )) is auto- 
matically self-consistent without the need to first solve 
a self-consistent integral equation, as is the case in the 
direct evaluation of (see Eq.(|lO|)). 

In practice, the direct evaluation of the time inte- 
gral in ( ^l|) is inconvenient because x{t) oscillates with 
undiminishing amplitude at large t. We have therefore 
noted that if we define a(u) = J °° X zz (t)e~ ut dt (i.e. 
a(u) — x zz (iu)), the VdW coefficient, Cq, can be written 
as 



C 6 = - 



du a A {u) a B (u). 



(25) 



For helium x(t) was calculated up to t — 15 atomic 
units (AU), which allows accurate calculation of a(u) for 
u > uo — 0.4. In the interval < u < uq, we repre- 
sented a(u) by the expression a + b/(l + cu 2 ), and fitted 
a, 6, c to a(u) and its first two derivatives at u — u a . 
(We checked that the results are insensitive to the ex- 
act choice of uq or to thw choice of the extrapolating 
function). Fig. 1 shows our a(u) for He. The correct 
asymoptotic form, a(u) — > 2/u 2 (the f-sum rule), is au- 
tomatically obeyed. The completeness sum rule, requires 
J °° a(u)du = 2n < o |z 2 |0o >— 2.50. Our a(u) gives 
2.33. An independent check on our a(u) is the static sus- 
ceptibility a(0). The best theoretical value is 1.383241 
Iffal, while we find 1.38. 
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Fig. 1. The imaginary-frequency susceptibility a(u) for Helium 
(solid line - direct evaluation, broken line - extrapolation). 

Our results for the He-He VdW constant is C§ = 1.45, 
almost identical to the best theoretical value Q 1.458. 
For the H-He system we find G§ — 2.81 compared to the 
best theoretical value ]17j of 2.817. 

We feel cautious about the significance of the high ac- 
curacy of our results for the He-He and the H-He systems 
in view of the fact that our calculated a(u) leads to a 7% 
error in the completeness sum rule. At the same time our 
results demonstrate the soundness and feasibility of our 
approach. We are optimistic that our approach will not 
only give asymptotic van der Waals coefficients, but the 
entire nuclear potential energy function e(R), including 
polarization energies. 

We found that the results are rather sensitive to the 
choice of a good KS potential for the unperturbed ground 
state. Repeating the calculation by replacing the exact 
V xc by V xc in the local density approximation, the result 
for Cq of the He-He system was 1.85, 28% too high. This 
is qualitatively similar to the experience of Petersilka et 
al. [O] with calculations of excited-state energies. 

We are indebted to C. Umrigar for providing us with 
the exact KS and the LDA KS data for Helium. We also 
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